Calcula la accesibilidad potencial que tienen tres municipios de una provincia de tu elección a mano de obra del entorno. Para ello, utiliza el indicador de potencial económico. La masa (M) será la población residente en todos los municipios de la provincia, puesto que el área de influencia no se limita a los tres municipios elegidos, y el coste (C) será la distancia que separa cada uno de los tres municipios de todos los demás municipios de la provincia calculada por la red de transporte.
El objetivo de este ejercicio es familiarizarse con la relación entre usos del suelo y transporte a través de un indicador de accesibilidad agregado. Para ello, se plantea calcular el potencial que tienen tres localizaciones para acceder a mano de obra. Por ejemplo, una gran industria puede tener como factor de localización una elevada masa de población residente lo más cerca posible. Desde el punto de vista técnico, el objetivo es dotar al estudiante de la habilidad de realizar un flujo de trabajo completo, desde la adquisición de datos, tratamiento, almacenamiento en una base de datos espacial, análisis y evaluación de resultados.
El desarrollo de este caso práctico consta de las siguientes tareas:
Se han utilizado los siguientes paquetes y versiones.
plic <- data.frame(installed.packages())
plic <- subset(plic, plic$Package%in% c('data.table', 'sp', 'rgdal', 'sf', 'leaflet','RColorBrewer','RPostgreSQL'))
plic <- plic[, c("Package", "Version")]
plic| Package | Version | |
|---|---|---|
| data.table | data.table | 1.14.8 |
| leaflet | leaflet | 2.2.0 |
| RColorBrewer | RColorBrewer | 1.1-3 |
| rgdal | rgdal | 1.6-7 |
| RPostgreSQL | RPostgreSQL | 0.7-5 |
| sf | sf | 1.0-14 |
| sp | sp | 2.0-0 |
## Loading required package: DBI## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE library(ggplot2) # gráficos
library(gridExtra) # capacidades extra ggplot
library(ggpmisc) # extensiones para ggplot2## Loading required package: ggpp##
## Attaching package: 'ggpp'## The following object is masked from 'package:ggplot2':
##
## annotate## Registered S3 method overwritten by 'ggpmisc':
## method from
## as.character.polynomial polynom## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal) library(leaflet.extras2) # capacidades adicionales
library(RColorBrewer) # paletas de colores
library(DT) # tablas dinámicas con datatable()
library(formattable) # tablas con estilo con formattable()
library(dplyr) # flujo de trabajo mejorado##
## Attaching package: 'dplyr'## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union library(sp) # para tratar datos espaciales como data.frames espaciales
library(rgdal) # para importar datos espaciales con sp## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/soho_/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:\Program Files\PostgreSQL\11\share\contrib\postgis-3.0\proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.Antes de comenzar, seguimos las instrucciones de los siguientes archivos para instalar PostgreSQL con las extensiones PostGIS y pgRouting, crear una base de datos, instalar y configurar el programa osm2po, y descargar la red de carreteras de España de OpenStreetMap, en nuestro caso de Cantabria:
01_Instalar_PostgreSQL_con_PostGIS.pdf 02_Crear_base_datos_pgAdmin.pdf 03_Descargar_carreteras_osm.pdf
En primer lugar, es necesario descargar e instalar PostgreSQL, y crear las extensiones PostGIS y pgRouting. Lo hacemos siguiendo las siguientes instrucciones: (instrucciones del archivo 01_Instalar_PostgreSQL_con_PostGIS.pdf)
Entramos en el siguiente enlace y descargamos la versión 11.21: https://www.enterprisedb.com/downloads/postgres-postgresql-downloads
Una vez descargado, hacemos doble click en el fichero, y esperamos unos segundos a que aparezca la pantalla de instalación para pinchar en siguiente.
Elegimos la ruta de instalación que aparece por defecto y seleccionamos todos los componentes por defecto.
Mantenemos el directorio de datos por defecto, hacemos click en “Siguiente” en las siguientes pantallas y esperamos mientras PostgreSQL se instala en nuestro ordenador (tarda unos minutos).
Escribimos una contraseña para el usuario “postgres”, y seleccionamos un puerto (por defecto el 5432).
Elegimos la configuración regional (por defecto) y le damos a “siguiente” al Resumen de la instalación.
Termina la instalación y lanzamos el Stack Builder.
Comprobamos que tenemos acceso a internet y elegimos la instalación sobre la que vamos a instalar las herramientas adicionales (PostgreSQL 11.21 en el puerto que hemos elegido 5432).
Desplegamos el menú “Spatial Extensions” y elegimos el Bundle de postGIS v.3.0.3.
Seleccionamos la ubicación para la descarga del archivo y comienza la instalación. Aceptamos los acuerdos de la licencia. Instalamos la extensión PostGIS y creamos una base de datos espacial.
Elegimos la carpeta de instalación, introducimos usuario “postgres”, contraseña y puerto de PostgreSQL “5432”.
Elegimos el nombre de la base de datos espacial. Es la base de datos que usaremos como plantilla para crear bases de datos espaciales, no debemos eliminarla: “postgis_30_sample”.
Aceptamos la instalación de GDAL_DATA para gestionar archivos raster.
Se requiere el reinicio de postgreSQL a continuación. Aceptamos el parámetro out de db rasters y queda la Instalación completada - Close.
Termina el Stack Builder, le damos a (Finish).
Abrimos pgAdmin (Inicio – pgAdmin 4) para comprobar que se han instalado las extensiones pgRouting y PostGIS.
Añadimos la ruta a psql.exe a las variables de entorno del sistema: En Windows, vamos al Panel de Control – Sistema – Configuración avanzada del sistema – Variables de entorno En el panel superior, seleccionamos PATH y pinchamos en Editar. En la nueva ventana, pinchamos en “Nueva” y pegamos la siguiente línea en la celda que se activa: C:Files\11 Pinchamos en “Aceptar” hasta que se cierren todas las ventanas y regresemos al Panel de Control.
En segundo lugar, crearemos una base de datos desde pgAdmin4 sobre la que vamos a trabajar, y la llamaremos “examen M12transportes”.
Una vez abierto el programa pgAdmin, desplegamos el elemento Servers y PostgreSQL 11 (puede pedir de nuevo la contraseña de la base de datos) y hacemos click con el botón derecho en Databases – Create – Database (asegurándonos de que la base de datos “postgis_30_sample” está desconectada (icono gris con aspa roja).
En la siguiente pantalla, damos un nombre a la base de datos, “examenM12”.
Pasa a la pestaña “Definition” y en “Template” elegimos la plantilla para la base de datos: “postgis_30_sample”. Se recomienda elegir el encoding UTF-8.
Por último, guardamos la base de datos haciendo click en “Save”.
Vamos a explicar las instrucciones para descargar la red de carreteras de “Cantabria” de Open Street Maps (OSM) y preparar la configuración para subirlo a una base de datos local PostGIS.
Descargamos el fichero con la red de carreteras de España. En el siguiente enlace están listados los servidores para la descarga de datos de OSM, con indicación de la cobertura espacial y la frecuencia de actualización: https://download.geofabrik.de/europe/spain.html
Descargamos el primer enlace de Commonly Used Formats, el fichero spain-latest.osm.pbf. En la dirección: “C:__4ºTRIMESTRE_TRANSPORTES_M12_osm-latest.osm.pbf”
Descargamos el programa “osm2po”. Usamos el programa libre osm2po para generar una sentencia SQL que permita subir la información del fichero de carreteras a la base de datos PostGIS. -Entramos en http://osm2po.de, usamos la última versión 5.5.5. -Guardamos el fichero, en una nueva carpeta “software” dentro del directorio EXAMEN_M12: “C:__4ºTRIMESTRE_TRANSPORTES_M122po-5.5.5”
Instalación de osm2po: Una vez descargado el fichero osm2po-5.1.0.zip, navegamos hasta la carpeta que lo contiene, descomprimimos en una carpeta con su mismo nombre y entramos en la nueva carpeta.
Accedemos a la configuración de la descarga: El siguiente paso es abrir el archivo “osm2po.config” para configurar la descarga de la zona y tipos de carreteras de Open Street Maps que nos convenga. Podemos editarlo con el Bloc de notas de Windows. Modificamos unas líneas para extraer la red de carreteras aptas para automóviles de España. Guardamos los cambios y cerramos el archivo.
Ahora ya tenemos el fichero, el software y la configuración necesaria.
Ya está todo listo para poder comenzar con el ejercicio, donde utilizaremos una red de transporte.
Ya podemos conectarnos desde R a nuestra base de datos PostgreSQL con funciones PostGIS y pgRouting.
# Setting connection parameters
host_db <- 'localhost'
port_db <- 5432
user_db <- 'postgres'
pass_db <- 'freskireski'
name_db <- 'examenM12'
max_db_conn <- 50
# Connecting to db
db_postgre_driver <- PostgreSQL(max.con = max_db_conn)
db_conn <- dbConnect(drv = db_postgre_driver,
user = user_db,
host = host_db,
port = port_db,
dbname = name_db,
pass = pass_db)Comprobamos que las extensiones han sido creadas
# https://rpostgres.r-dbi.org/reference/postgres-query.html
dbGetQuery(db_conn, "SELECT extname, extversion FROM pg_extension;")| extname | extversion |
|---|---|
| plpgsql | 1.0 |
| postgis | 3.0.3 |
| postgis_raster | 3.0.3 |
| pgrouting | 3.1.1 |
| postgis_topology | 3.0.3 |
| address_standardizer | 3.0.3 |
| postgis_sfcgal | 3.0.3 |
| pointcloud | 1.2.1 |
| pointcloud_postgis | 1.2.1 |
| ogr_fdw | 1.0 |
| fuzzystrmatch | 1.1 |
| postgis_tiger_geocoder | 3.0.3 |
Los archivos que se utilizan a continuación, están en “C:__4ºTRIMESTRE_TRANSPORTES_M12_data" Primero, establecemos el directorio de trabajo.
# directorio de trabajo.
wd_02 <- "C:/Users/soho_/Desktop/DATOS_4ºTRIMESTRE/M12_TRANSPORTES/EXAMEN_M12"
# Setting input_path
input_path_02 <- paste0(wd_02, "/input_data/")Vamos a importar un fichero shp con la geometría de los municipios de Cantabria, establecemos el sistema de coordenadas que indica la fuente de datos, y lo transformamos al sistema de coordenadas común con el resto de fuentes.
Fuente: Líneas límite municipales, Información geográfica de referencia, IGN. http://centrodedescargas.cnig.es/CentroDescargas/index.jsp.
GEOMETRÍA MUNICIPIOS shp
# Reading shp file as sf
fname_municipios <- paste0(input_path_02, "municipios_ign/lineas_limite/SIGLIM_Publico_INSPIRE/SHP_ETRS89/recintos_municipales_inspire_peninbal_etrs89/recintos_municipales_inspire_peninbal_etrs89.shp")
municipios_sf <- sf::st_read(fname_municipios,
type = 6, # Multipolygon
stringsAsFactors = FALSE,
int64_as_string = FALSE) # EPSG 4258## Reading layer `recintos_municipales_inspire_peninbal_etrs89' from data source
## `C:\Users\soho_\Desktop\DATOS_4ºTRIMESTRE\M12_TRANSPORTES\EXAMEN_M12\input_data\municipios_ign\lineas_limite\SIGLIM_Publico_INSPIRE\SHP_ETRS89\recintos_municipales_inspire_peninbal_etrs89\recintos_municipales_inspire_peninbal_etrs89.shp'
## using driver `ESRI Shapefile'
## converted into: MULTIPOLYGON
## Simple feature collection with 8129 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -9.301516 ymin: 35.17045 xmax: 4.327785 ymax: 43.79238
## Geodetic CRS: ETRS89# See https://en.wikipedia.org/wiki/Well-known_text#Well-known_binary for data types
head(municipios_sf)| INSPIREID | COUNTRY | NATLEV | NATLEVNAME | NATCODE | NAMEUNIT | CODNUT1 | CODNUT2 | CODNUT3 | geometry |
|---|---|---|---|---|---|---|---|---|---|
| ES.IGN.SIGLIM34091717035 | ES | https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder | Municipio | 34091717035 | Camós | ES5 | ES51 | ES512 | MULTIPOLYGON (((2.732126 42… |
| ES.IGN.SIGLIM34091717036 | ES | https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder | Municipio | 34091717036 | Campdevànol | ES5 | ES51 | ES512 | MULTIPOLYGON (((2.09486 42…. |
| ES.IGN.SIGLIM34091717037 | ES | https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder | Municipio | 34091717037 | Campelles | ES5 | ES51 | ES512 | MULTIPOLYGON (((2.097592 42… |
| ES.IGN.SIGLIM34091717038 | ES | https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder | Municipio | 34091717038 | Campllong | ES5 | ES51 | ES512 | MULTIPOLYGON (((2.816395 41… |
| ES.IGN.SIGLIM34091717039 | ES | https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder | Municipio | 34091717039 | Camprodon | ES5 | ES51 | ES512 | MULTIPOLYGON (((2.31534 42…. |
| ES.IGN.SIGLIM34091717040 | ES | https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder | Municipio | 34091717040 | Canet d’Adri | ES5 | ES51 | ES512 | MULTIPOLYGON (((2.677743 42… |
# Keeping Cantabria municipios (CODNUT3 = ES130)
# NUTs codes: https://es.wikipedia.org/wiki/NUTS_de_España
municipios_sf <- subset(municipios_sf, municipios_sf$CODNUT3 == "ES130")
# Setting CRS (ETRS89 long,lat [degree]) No UTM zone
municipios_sf <- sf::st_sf(municipios_sf,
crs = 4258)
# Transforming to EPSG:4326
municipios_sf <- sf::st_transform(municipios_sf,
crs = 4326)POBLACIÓN INE
Importamos un fichero csv con la población por municipios.
Fuente: Padrón municipal de habitantes del INE, https://www.ine.es/dynt3/inebase/index.htm?padre=525.
# Reading csv file as data.table
filename <- paste0(input_path_02, "poblacion_ine/ine_cantabria.csv")
poblacion_dt <- data.table::fread(filename, sep = ";", dec = ',', encoding = 'Latin-1')
head(poblacion_dt)| Municipios | Sexo | Periodo | Total |
|---|---|---|---|
| 39 Cantabria | Total | 2022 | 585.402 |
| 39 Cantabria | Total | 2021 | 584.507 |
| 39 Cantabria | Total | 2020 | 582.905 |
| 39 Cantabria | Total | 2019 | 581.078 |
| 39 Cantabria | Total | 2018 | 580.229 |
| 39 Cantabria | Total | 2017 | 580.295 |
## Classes 'data.table' and 'data.frame': 8343 obs. of 4 variables:
## $ Municipios: chr "39 Cantabria" "39 Cantabria" "39 Cantabria" "39 Cantabria" ...
## $ Sexo : chr "Total" "Total" "Total" "Total" ...
## $ Periodo : int 2022 2021 2020 2019 2018 2017 2016 2015 2014 2013 ...
## $ Total : chr "585.402" "584.507" "582.905" "581.078" ...
## - attr(*, ".internal.selfref")=<externalptr>| Municipios | Sexo | Periodo | Total |
|---|---|---|---|
| 39 Cantabria | Total | 2022 | 585.402 |
| 39 Cantabria | Hombres | 2022 | 283.718 |
| 39 Cantabria | Mujeres | 2022 | 301.684 |
| 39001 Alfoz de Lloredo | Total | 2022 | 2.466 |
| 39001 Alfoz de Lloredo | Hombres | 2022 | 1.284 |
| 39001 Alfoz de Lloredo | Mujeres | 2022 | 1.182 |
Transformamos la tabla para tener una columna por cada atributo de población (Total, Hombres, Mujeres) e ignorando el año (que es siempre el mismo 2022)
# Removing thousands separator and converting to integer type
poblacion_dt22 <- poblacion_dt22[, Total := as.integer(gsub('.', '', Total, fixed = TRUE))]
# Transposing the table
poblacion_dt22 <- dcast.data.table(poblacion_dt22, Municipios ~ Sexo, value.var = 'Total')
head(poblacion_dt22)| Municipios | Hombres | Mujeres | Total |
|---|---|---|---|
| 39 Cantabria | 283718 | 301684 | 585402 |
| 39001 Alfoz de Lloredo | 1284 | 1182 | 2466 |
| 39002 Ampuero | 2196 | 2188 | 4384 |
| 39003 Anievas | 147 | 120 | 267 |
| 39004 Arenas de Iguña | 870 | 864 | 1734 |
| 39005 Argoños | 945 | 903 | 1848 |
CREAMOS UNIÓN DE LAS DOS TABLAS MUNICIPIOS+POBLACIÓN
Tratamos los campos de código o nombre del municipio para poder añadir la información de población a los polígonos de los municipios.
Necesitamos los los 5 primeros dígitos del campo Municipios de la tabla con la población y los 5 últimos dígitos del campo NATCODE de la capa de municipios.
# Creating codigo_ine field
poblacion_dt22$codigo_ine <- substring(poblacion_dt22$Municipios, 1, 5)
municipios_sf$codigo_ine <- substring(municipios_sf$NATCODE, 7, 11)
# Adding poblacion_dt fields to municipios by codigo_ine
municipios_sf <- merge(municipios_sf,
poblacion_dt22,
by = 'codigo_ine',
all.x = TRUE)
# calculamos la superficie y densidad de habitantes
municipios_sf$superficie <- round(st_area(municipios_sf)/10^6,2) # m^2 a km^2
municipios_sf$densidad <- round(municipios_sf$Total / municipios_sf$superficie, 2)
# quitamos unidades (para ggplot es necesario)
municipios_sf$densidad <- units::drop_units(municipios_sf$densidad)
municipios_sf$superficie <- units::drop_units(municipios_sf$superficie)Comprobamos que tenemos nuestros datos bien descritos.
##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rows## Classes 'sf' and 'data.frame': 103 obs. of 17 variables:
## $ codigo_ine: chr "39001" "39002" "39003" "39004" ...
## $ INSPIREID : chr "ES.IGN.SIGLIM34063939001" "ES.IGN.SIGLIM34063939002" "ES.IGN.SIGLIM34063939003" "ES.IGN.SIGLIM34063939004" ...
## $ COUNTRY : chr "ES" "ES" "ES" "ES" ...
## $ NATLEV : chr "https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder" "https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder" "https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder" "https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/4thOrder" ...
## $ NATLEVNAME: chr "Municipio" "Municipio" "Municipio" "Municipio" ...
## $ NATCODE : chr "34063939001" "34063939002" "34063939003" "34063939004" ...
## $ NAMEUNIT : chr "Alfoz de Lloredo" "Ampuero" "Anievas" "Arenas de Iguña" ...
## $ CODNUT1 : chr "ES1" "ES1" "ES1" "ES1" ...
## $ CODNUT2 : chr "ES13" "ES13" "ES13" "ES13" ...
## $ CODNUT3 : chr "ES130" "ES130" "ES130" "ES130" ...
## $ Municipios: chr "39001 Alfoz de Lloredo" "39002 Ampuero" "39003 Anievas" "39004 Arenas de Iguña" ...
## $ Hombres : int 1284 2196 147 870 945 1092 271 8861 2316 335 ...
## $ Mujeres : int 1182 2188 120 864 903 1078 200 9292 2187 328 ...
## $ Total : int 2466 4384 267 1734 1848 2170 471 18153 4503 663 ...
## $ geometry :sfc_MULTIPOLYGON of length 103; first list element: List of 1
## ..$ :List of 1
## .. ..$ : num [1:704, 1:2] -4.27 -4.26 -4.26 -4.26 -4.25 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
## $ superficie: num 46.69 32.41 21.62 86.04 5.56 ...
## $ densidad : num 52.8 135.3 12.3 20.1 332.4 ...
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:16] "codigo_ine" "INSPIREID" "COUNTRY" "NATLEV" ...# tabla
tabla<-municipios_sf %>%as.data.frame %>% select (codigo_ine, CODNUT3, Municipios, Hombres, Mujeres, Total, superficie, densidad )
tabla%>%datatable()library(mapview)
mapviewOptions(fgb = FALSE) # needed when creating web pages
pink2 = colorRampPalette(c('lightgreen', 'darkgreen'))
mapview(municipios_sf[c("Municipios","Total","densidad")],
# alpha.regions = 1,
map.types = c("OpenStreetMap.DE", "CartoDB.Positron","Esri.WorldShadedRelief"),
layer.name = "Densidad",
zcol = "densidad",
col.regions = pink2, at = (seq(0, 400, 50)),
fgb = FALSE)Vamos a usar el programa osm2po para generar un fichero sql con la sentencia para subir la red de carreteras que hemos descargado siguiendo las instrucciones del archivo 03_Descargar_carreteras_osm.pdf a nuestra base de datos PostGIS.
Lo vamos a hacer usando la función system para ejecutar osm2po sin salir de R.
Cuando aparezca el siguiente mensaje en la consola puedes parar la ejecución presionando el icono Stop de la consola.
INFO WebDir is C:\Users\soho_\Desktop\DATOS_4ºTRIMESTRE\M12_TRANSPORTES\EXAMEN_M12\software\osm2po-5.5.5\osm2po-web
INFO https://localhost:8889/Osm2poService/ STARTED - 221022-10:51.16360 - 3.759M
INFO http://localhost:8888/Osm2poService/ STARTED - 221022-10:51.16362 - 3.758M
WARN Log statistics: WARNs:2
INFO All Services STARTED. Press Ctrl-C to stop them.En ese momento y se habrá generado el fichero sp2_2po_4pgr.sql dentro de la carpeta *“C:__4ºTRIMESTRE_TRANSPORTES_M122po-5.5.5“*.
Comprobamos si la topología es correcta con la función st_is_valid.
# Checking topology
topo_check <- sf::st_is_valid(municipios_sf)
# Percentage of elements with valid topology
topo_perc_ok <- length(subset(topo_check, topo_check == TRUE)) * 100 / nrow(municipios_sf)
# Number of elements with invalid topology
topo_number_ko <- length(subset(topo_check, topo_check == FALSE))
# id of elements with invalid topology
topo_elements_not_ok <- subset(topo_check, topo_check == FALSE)
# Returns
topo_results <- c("topo_perc_ok" = topo_perc_ok,
"topo_number_ko" = topo_number_ko,
"topo_elements_not_ok" = topo_elements_not_ok)
topo_results## topo_perc_ok topo_number_ko
## 100 0Comprobamos que tiene una geometría válida.
SUBIMOS A POSTGIS
Subimos la capa a PostGIS. Para que no haya problemas con los nombres de las columnas cuando hagamos las consultas, antes de subir la capa convertimos todos los nombres a minúsculas.
# Setting column names to lower case
colnames(municipios_sf) <- tolower(colnames(municipios_sf))
# Writing municipios in PostGIS table
sf::st_write(dsn = db_conn,
obj = municipios_sf,
layer = "municipios_cantabria",
layer_options = c("OVERWRITE=yes"))
# Setting SRID
srid_query <- "
SELECT UpdateGeometrySRID('municipios_cantabria','geometry',4326);
"
dbExecute(db_conn, srid_query)Creamos una capa con los centroides del municipio
(en realidad punto similar al centroide que cae dentro del polígono municipal correspondiente).
Esta operación la realizamos directamente en la base de datos, utilizando la función ST_PointOnSurface sobre el campo geometry, que es el que almacena la información espacial.
# Creating a point layer
# ST_PointOnSurface Returns a POINT guaranteed to lie on the surface
# https://postgis.net/docs/ST_PointOnSurface.html
centroids_query <- paste0("
DROP TABLE IF EXISTS municipios_cantabria_point;
CREATE TABLE municipios_cantabria_point AS
SELECT codigo_ine,
inspireid,
country,
nameunit,
codnut1,
codnut2,
codnut3,
total,
hombres,
mujeres,
ST_PointOnSurface(geometry) as geometry
FROM municipios_cantabria;
")
dbExecute(db_conn, centroids_query) ## [1] 103# Setting SRID
srid_query <- "
SELECT UpdateGeometrySRID('municipios_cantabria','geometry',4326);
"
dbExecute(db_conn, srid_query)## [1] -1Subimos la red de carreteras a PostGIS para poder crear una red de transporte a partir del fichero sp2_2po_4pgr.sql generado anteriormente. Para ello ejecutamos psql desde R usando la función shell.
Para que este comando funcione, es necesario que esté definida la ruta (PATH) de psql en "Variable de entornos". En nuestro caso: C:/Program Files/PostgreSQL/11/binPero tenemos que ejecutarlo en el shell de windows porque no se ejecuta:
Uploading sql file to PostGISsystem(‘psql -U postgres -d examenM12 –host localhost –port 5432 -q -f “C:/Users/soho_/Desktop/DATOS_4ºTRIMESTRE/M12_TRANSPORTES/EXAMEN_M12/software/osm2po-5.5.5/sp2/sp2_2po_4pgr.sql”’)
addgeometrycolumn
-------------------------------------------------------------
public.sp2_2po_4pgr.geom_way SRID:4326 TYPE:LINESTRING DIMS:2
(1 fila)Para reducir el tiempo de computación de las distancias a través de la red de transporte, vamos a crear una tabla con la red de carreteras de Cantabria y alrededores. Para ello, vamos a usar las funciones ST_Extent (para calcular el rectángulo que encuadra los municipios de Cantabria) y ST_Buffer (para añadir unos km más alrededor del rectángulo y asegurar la conexión por carretera).
# Creating bbox_cantabria table
query_bbox <- "
DROP TABLE IF EXISTS bbox_cantabria;
CREATE TABLE bbox_cantabria AS
SELECT ST_SetSRID(ST_Extent(ST_Buffer(geometry, 0.1)),4326) as bextent
FROM municipios_cantabria;
"
dbExecute(db_conn, query_bbox)## [1] 1# Creating carreteras_cantabria table
query_carreteras <- "
DROP TABLE IF EXISTS carreteras_cantabria;
CREATE TABLE carreteras_cantabria AS (
SELECT
id,
osm_id,
osm_name,
source,
target,
km,
kmh,
cost,
reverse_cost,
x1,
y1,
x2,
y2,
geom_way
FROM sp2_2po_4pgr
WHERE ST_INTERSECTS(geom_way,
(SELECT bextent
FROM bbox_cantabria))
);
"
dbExecute(db_conn, query_carreteras)## [1] 174247El siguiente paso es crear la topología de la red, para lo que necesitamos crear la tabla de vértices únicos (a partir de los campos source y target).
# If the road network table already has source and target fields
# with unique vertex identifiers assigned to the corresponding edges,
# then we only need to create the vertices table, which should be fast (less than 1 minute)
# Creating vertices table
# https://docs.pgrouting.org/2.6/en/pgr_createVerticesTable.html
create_vertices_table_query <- "
SELECT pgr_createVerticesTable('carreteras_cantabria',
the_geom:='geom_way',
source:='source',
target:='target'
);"
dbGetQuery(db_conn, create_vertices_table_query)| pgr_createverticestable |
|---|
| OK |
# However, if we do not have these fields properly filled, we need to add network topology
# which is a much resource-consuming task
# Info: https://blog.daftcode.pl/find-your-way-with-the-power-of-postgis-pgrouting-66d620ef201bAntes de seguir, es conveniente analizar la topología creada.
Seleccionamos el nodo de de la red de carreteras que es más cercano al centroide de cada municipio. Este nodo será el verdadero punto de origen y destino en los cálculos.
nearest_vertex_query <- "
DROP TABLE IF EXISTS municipio_nearest_vertex;
CREATE TABLE municipio_nearest_vertex AS
SELECT munis.codigo_ine,
distance,
ver_id
FROM municipios_cantabria_point munis
JOIN LATERAL (
SELECT
m.codigo_ine,
st_distance(ST_Transform(m.geometry, 25830),ST_Transform(v.the_geom, 25830)) AS distance,
v.id as ver_id
FROM municipios_cantabria_point m, carreteras_cantabria_vertices_pgr v
WHERE codigo_ine = munis.codigo_ine
ORDER BY m.geometry <-> v.the_geom
limit 1
) p on true
ORDER BY munis.codigo_ine;
"
dbExecute(db_conn, nearest_vertex_query)## [1] 103Calculamos la ruta más corta por la red de transporte entre “Santander”(código ine:39075) y “San Vicente de la Barquera”(código ine:39080) utilizando el algoritmo de Dijkstra.
# Shortest route between two points
# https://docs.pgrouting.org/latest/en/pgr_dijkstra.html
# pgr_dijkstra(edges_sql, start_vid, end_vid [, directed])
shortest_route_query <- "SELECT * FROM pgr_dijkstra(
'SELECT id,
source,
target,
km * 1000 as cost,
km * 1000 as reverse_cost
FROM carreteras_cantabria',
(SELECT ver_id
FROM municipio_nearest_vertex
WHERE codigo_ine::int = 39075),
(SELECT ver_id
FROM municipio_nearest_vertex
WHERE codigo_ine::int = 39080)
);"
shortest_route <- data.table(dbGetQuery(db_conn, shortest_route_query))
shortest_route%>%datatable()Seleccionamos la distancia total:
# Selecting the value of agg_cost column on the last row (where column edge value is -1), rounding to 3 decimal digits.
distancia_total <- round(shortest_route[edge == -1, .(agg_cost)], 3)
print(paste0("Distancia total de la ruta: ", distancia_total, " metros."))## [1] "Distancia total de la ruta: 56138.75 metros."Vemos que hay 56 kilómetros aproximadamente entre los municipios “Santander”(código ine:39075) y “San Vicente de la Barquera”(código ine:39080).
Calculamos la matriz de distancias entre todos los municipios de Cantabria (todos contra todos) utilizando la función pgr_dijkstraCostMatrix.
# Vértices más cercanos
vertex_query <- "
SELECT codigo_ine,
ver_id
FROM municipio_nearest_vertex
WHERE codigo_ine IN (SELECT codigo_ine
FROM municipios_cantabria_point
WHERE total IS NOT NULL
ORDER BY total DESC
);
"
vertex_selection <- data.table(dbGetQuery(db_conn, vertex_query))
# Matriz de costes de Dijkstra
dijkstra_route_query <- "
SELECT * FROM pgr_dijkstraCostMatrix(
'SELECT id,
source,
target,
km * 1000 as cost,
km * 1000 as reverse_cost
FROM carreteras_cantabria',
(SELECT array_agg(id)
FROM carreteras_cantabria_vertices_pgr
WHERE id IN
(SELECT ver_id
FROM municipio_nearest_vertex
WHERE codigo_ine IN (SELECT codigo_ine
FROM municipios_cantabria_point
WHERE total IS NOT NULL
ORDER BY total DESC
)
)
),
directed := false
);"
matriz_origen_destino <- data.table(dbGetQuery(db_conn, dijkstra_route_query))
# Añadimos el codigo_ine inicial al resultado
matriz_origen_destino <- merge(x = matriz_origen_destino,
y = vertex_selection,
by.x = 'start_vid',
by.y = 'ver_id')
colnames(matriz_origen_destino) <- c('star_vid', 'end_vid', 'agg_cost', 'codigo_origen')
# Añadimos el codigo_ine final al resultado
matriz_origen_destino <- merge(x = matriz_origen_destino,
y = vertex_selection,
by.x = 'end_vid',
by.y = 'ver_id')
colnames(matriz_origen_destino) <- c('end_vid', 'start_vid', 'agg_cost', 'codigo_origen', 'codigo_destino')
# Visualizamos los resultados
matriz_origen_destino%>%datatable()la tabla indica los vértices de entrada y salida de cada tramo, el costo acumulado (distancia entre municipios) y los códigos_ine de los municipios.
Se nos plantea un problema en el que tenemos que elegir tres municipios, así que creamos una tabla con sus códigos INE:
municipios_examen <- data.frame(municipios = c("Santander",
"San Vicente de la Barquera",
"Torrelavega"),
codigo_ine = c("39075",
"39080",
"39087"))
municipios_examen %>% formattable()| municipios | codigo_ine |
|---|---|
| Santander | 39075 |
| San Vicente de la Barquera | 39080 |
| Torrelavega | 39087 |
Se trataría ahora de filtrar por los códigos ine de los municipios elegidos: Santander (39075), San Vicente de la Barquera (39080) y Torrelavega (39087), de la matriz origen_destino.
MODexamen<-matriz_origen_destino %>% filter(codigo_origen %in% municipios_examen$codigo_ine) %>% datatable()
MODexamenTenemos un total de 102 municipios. Cada uno de los tres municipios tiene 101 conexiones. Esto hace un total de 3x101 = 303 conexiones que son las que se muestran en la tabla.
El modelo de potencial económico es sin duda la medida gravitatoria más utilizada en los estudios de accesibilidad (ver por ejemplo Hansen, 1959 o Dundon-Smith y Gibb, 1994). Este modelo se basa en que la accesibilidad de un nodo dado en una red es proporcional a la interacción espacial entre ese nodo y los demás nodos de la red. El grado de oportunidades entre dos áreas está positivamente relacionado con la capacidad de atracción de las áreas y negativamente relacionado con la impedancia entre las áreas (Linneker y Spence, 1992).
Por ejemplo, la accesibilidad de una ciudad es proporcional al tamaño de las ciudades que la rodean (en población, centros de empleo, etc.) e inversamente proporcional a la distancia (o tiempo o coste de viaje) entre esa ciudad y las que la rodean. Así, el potencial económico del nodo i, Pi, se obtiene según la siguiente ecuación:
Pi=∑Mj/exp(Cij)α
donde:
Cij es la distancia entre los nodos (centroides) i y j , Mj es la
masa (población) de los municipios
el parámetro α,es un parámetro que suele oscilar entre 1.0 y 2.5 (cuanto
mayor sea el valor de α mayor la penalización por distancia). En la
mayor parte de los estudios empíricos toma el valor α= 1.
El potencial económico tiene unidades y depende de si medimos la distancia en metros o kilómetros tendremos una magnitud mayor o menor. Aquí decido dejar la distancia en metros.
Los datos que necesitamos, por tanto, para este análisis es distancia y población. Esto es: una matriz de origen destino que ya tenemos calculada del paso anterior donde, con la función merge podemos agregar la población que tenemos en la tabla “poblacion_dt” gracias a la columna del código INE común en las dos.
Matriz de Coste-Población para Santander (39075), San Vicente de la Barquera (39080) y Torrelavega (39087):
# matriz de coste-población
matriz_potencial <- merge(x = matriz_origen_destino%>% filter(codigo_origen %in% municipios_examen$codigo_ine),
y = poblacion_dt22,
by.x = 'codigo_destino',
by.y = 'codigo_ine') %>%
select(codigo_origen, codigo_destino, agg_cost, Municipios, Total)
matriz_potencial <- matriz_potencial[, c(1,2,3,5)]
# Ordenamos por el código ine del municipio origen
matriz_potencial=matriz_potencial[order(codigo_origen)]
matriz_potencial %>% datatable()Definimos una función que nos calcula el potencial en función de la matriz que tenemos, el código ine del municipio y el valor de α deseado:
# Función para calcular el potencial económico de un municipio
# identificado por su codigo_ine
# dada la matriz de coste-poblacion de una región
potencial_economico <- function(matriz_datos, ine_municipio, alpha = 1){
potencial_economico <- matriz_datos %>%
filter(codigo_origen == ine_municipio) %>%
mutate(potencial_unitario = Total / agg_cost^alpha)
potencial_economico <- round(sum(potencial_economico$potencial_unitario), 3)
}Podemos utilizarla para calcular el potencial económico de cualquiera de los municipios elegidos para cualquier alfa. Tomemos por ejemplo el valor de alfa = 1. Otra de las mejoras, sería modificar el parámetro α: El parámetro α impone una penalización a las distancias. Usando la fórmula anterior se han analizado tres valores de α: 1, 1.5 y 1.7. Los resultados obtenidos son los siguientes:
POTENCIAL CON ALFA=1
potencial_municipios <- sapply(municipios_examen$codigo_ine,
FUN = potencial_economico,
matriz_datos = matriz_potencial,
alpha = 1,
USE.NAMES = TRUE)
potencial_municipios <- as.data.frame(potencial_municipios)
potencial_municipios$codigo_ine <- rownames(potencial_municipios)
rownames(potencial_municipios) <- NULL
potencial_municipios <- merge(x = potencial_municipios,
y = poblacion_dt22[,c("codigo_ine","Municipios")],
by.x = "codigo_ine",
by.y = "codigo_ine")
colnames(potencial_municipios) <- c("codigo_ine", "potencial", "municipios")
# reordenamos las columnas
potencial_municipios <-potencial_municipios[,c(3,1,2)] %>% arrange(desc(potencial))
potencial_municipios %>% formattable()| municipios | codigo_ine | potencial |
|---|---|---|
| 39087 Torrelavega | 39087 | 23.660 |
| 39075 Santander | 39075 | 19.059 |
| 39080 San Vicente de la Barquera | 39080 | 12.161 |
POTENCIAL CON ALFA=1.5
potencial_municipios2 <- sapply(municipios_examen$codigo_ine,
FUN = potencial_economico,
matriz_datos = matriz_potencial,
alpha = 1.5,
USE.NAMES = TRUE)
potencial_municipios2 <- as.data.frame(potencial_municipios2)
potencial_municipios2$codigo_ine <- rownames(potencial_municipios2)
rownames(potencial_municipios2) <- NULL
potencial_municipios2 <- merge(x = potencial_municipios2,
y = poblacion_dt22[,c("codigo_ine","Municipios")],
by.x = "codigo_ine",
by.y = "codigo_ine")
colnames(potencial_municipios2) <- c("codigo_ine", "potencial", "municipios")
# reordenamos las columnas
potencial_municipios2 <-potencial_municipios2[,c(3,1,2)] %>% arrange(desc(potencial))
potencial_municipios2 %>% formattable()| municipios | codigo_ine | potencial |
|---|---|---|
| 39087 Torrelavega | 39087 | 0.178 |
| 39075 Santander | 39075 | 0.153 |
| 39080 San Vicente de la Barquera | 39080 | 0.062 |
Sigue teniendo mejor potencial Torrelavega.
POTENCIAL ALFA=1.7
potencial_municipios3 <- sapply(municipios_examen$codigo_ine,
FUN = potencial_economico,
matriz_datos = matriz_potencial,
alpha = 1.7,
USE.NAMES = TRUE)
potencial_municipios3 <- as.data.frame(potencial_municipios3)
potencial_municipios3$codigo_ine <- rownames(potencial_municipios3)
rownames(potencial_municipios3) <- NULL
potencial_municipios3 <- merge(x = potencial_municipios3,
y = poblacion_dt22[,c("codigo_ine","Municipios")],
by.x = "codigo_ine",
by.y = "codigo_ine")
colnames(potencial_municipios3) <- c("codigo_ine", "potencial", "municipios")
# reordenamos las columnas
potencial_municipios3 <-potencial_municipios3[,c(3,1,2)] %>% arrange(desc(potencial))
potencial_municipios3 %>% formattable()| municipios | codigo_ine | potencial |
|---|---|---|
| 39087 Torrelavega | 39087 | 0.026 |
| 39075 Santander | 39075 | 0.023 |
| 39080 San Vicente de la Barquera | 39080 | 0.008 |
Se siguen manteniendo los mismos potenciales relativos entre los tres municipios. Siendo el mayor potencial el de Torrelavega.
# preparamos la tabla
alphas <- c(1, 1.5, 1.7)
resultados <- rbind(potencial_municipios, potencial_municipios2, potencial_municipios3)
resultados <- cbind(alpha = alphas, resultados)
rownames(resultados) <- NULL
resultados<-data.table(resultados)
formattable(resultados)| alpha | municipios | codigo_ine | potencial |
|---|---|---|---|
| 1.0 | 39087 Torrelavega | 39087 | 23.660 |
| 1.5 | 39075 Santander | 39075 | 19.059 |
| 1.7 | 39080 San Vicente de la Barquera | 39080 | 12.161 |
| 1.0 | 39087 Torrelavega | 39087 | 0.178 |
| 1.5 | 39075 Santander | 39075 | 0.153 |
| 1.7 | 39080 San Vicente de la Barquera | 39080 | 0.062 |
| 1.0 | 39087 Torrelavega | 39087 | 0.026 |
| 1.5 | 39075 Santander | 39075 | 0.023 |
| 1.7 | 39080 San Vicente de la Barquera | 39080 | 0.008 |
Se comprueba que el municipio con más potencial para una actividad con alta demanda de mano de obra, sería Torrelavega seguido de Santander. Y en último lugar, sería San Vicente de la Barquera el municipio elegido para crear la actividad con alta demanda de empleo.
Veámoslo geográficamente:
potencial_sf <- municipios_sf[,c("codigo_ine", "geometry", "Municipios")]
potencial_sf <- merge(potencial_sf,
potencial_municipios[,c("codigo_ine", "potencial")],
by = 'codigo_ine',
all.x = TRUE)
mapview(potencial_sf,
alpha.regions = 1,
map.types = 'CartoDB.Positron',
layer.name = "potencial",
zcol = "potencial",
fgb = FALSE)De la misma manera es posible calcular el potencial para todos los municipios. Usando el valor de alfa= 1, tendríamos estos valores:
# matriz de coste-población-todos municipios
matriz_potencial_todos <- merge(x = matriz_origen_destino,
y = poblacion_dt22,
by.x = 'codigo_destino',
by.y = 'codigo_ine') %>%
select(codigo_origen, codigo_destino, agg_cost, Municipios, Total)
matriz_potencial_todos <- matriz_potencial_todos[, c(1,2,3,5)]
# Ordenamos por el código ine del municipio origen
matriz_potencial_todos=matriz_potencial_todos[order(codigo_origen)]
# Función para calcular el potencial económico de un municipio
# identificado por su codigo_ine
# dada la matriz de coste-poblacion de una región
potencial_economico <- function(matriz_datos, ine_municipio, alpha = 1){
potencial_economico <- matriz_datos %>%
filter(codigo_origen == ine_municipio) %>%
mutate(potencial_unitario = Total / agg_cost^alpha)
potencial_economico <- round(sum(potencial_economico$potencial_unitario), 3)
}
potencial_municipios_todos <- sapply(municipios_sf$codigo_ine,
FUN = potencial_economico,
matriz_datos = matriz_potencial_todos,
alpha = 1,
USE.NAMES = TRUE)
potencial_municipios_todos<- as.data.frame(potencial_municipios_todos)
potencial_municipios_todos$codigo_ine <- rownames(potencial_municipios_todos)
rownames(potencial_municipios_todos) <- NULL
potencial_municipios_todos <- merge(x = potencial_municipios_todos,
y = poblacion_dt22[,c("codigo_ine","Municipios")],
by.x = "codigo_ine",
by.y = "codigo_ine")
colnames(potencial_municipios_todos) <- c("codigo_ine", "potencial", "municipios")
# reordenamos las columnas
potencial_municipios_todos <- potencial_municipios_todos[,c(3,1,2)] %>% arrange(desc(potencial))
potencial_municipios_todos%>%datatable()CONCLUSIÓN Los tres municipios donde sería óptimo crear una actividad que demande mano de obra, serían por orden de prioridad:
Santa Cruz de Bezana, El Astillero y Camargo.
¿Qué datos aportarían mayor nivel de detalle espacial/temporal? Entre las mejoras que podríamos incorporar sería obtener la población que está desempleada y en edad de trabajar, para completar y garantizar que la mano de obra necesaria se pueda obtener. Esta información se encuentra en el INE.
¿Qué sucedería si se incrementa o disminuye el parámetro de α de decaimiento con la distancia? Y la otra mejora ya la hemos introducido en el estudio, al probar con diferentes valores del parámetro “alfa” (1,1.5,1.7).